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Abstract. The approach to the singularity in Gowdy spacetimes consists of velocity 
term dominated behavior, except at a set of isolated points. At and near these points, 
spiky features grow. This paper reviews what is known about these spikes. 
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1. Introduction 

There have been several investigations, both analytical and numerical, of the approach 
to the singularity in inhomogeneous cosmologies. P^-JI] The most extensively studied 
class of these spacetimes is the class of Gowdy spacetimes [IH] on T 3 x R. Theorems 
due to Isenberg and MoncriefJT] showed that the polarized Gowdy spacetimes are 
asymptotically velocity term dominated (AVTD). It was then natural to ask whether this 
is also true of the more general unpolarized Gowdy spacetimes. This was a more difficult 
question since the relevant equations are linear in the polarized case but nonlinear in the 
unpolarized case. To address this question, Berger and Moncrief performed numerical 
simulations of the approach to the singularity in Gowdy spacetimes. They obtained the 
following curious result: the behavior is AVTD except at a set of isolated points. At and 
near these exceptional points, features appear that grow ever steeper and narrower as 
the singularity is approached. Berger and Moncrief referred to these features as "spiky 
features," though for brevity we will simply call them "spikes." These spikes form a 
"fine structure" on the simpler and more prevalent behavior of the rest of the spacetime. 

The purpose of this paper is to review what is known about these spikes. Section 
2 presents the basic equations governing the Gowdy spacetimes on T 3 x R. The next 
three sections concern a variety of methods and results, with numerical methods and 
results presented in section 3, analytical approximations in section 4 and mathematical 
results in section 5. Conclusions are given in section 6. 

2. Equations 

The Gowdy metric on T 3 x R takes the form 

ds 2 = e x/2 r 1/2 (-dt 2 + dx 2 ) + t[e p (dy + Qdz) 2 + e~ p dz 2 ] (1) 

where P, Q and A are functions of t and x. The T 3 spatial topology is imposed by having 
< x, y, z < 2tc and having P, Q and A be periodic functions of x. The vacuum Einstein 
field equations split into "evolution" equations for P and Q 

P,tt + r 1 ^ - P, xx + e 2P (Q 2 x - Q 2 t ) = (2) 

Q,tt + t^Q,t - Q, xx + 2(P >t Q jt - P tX Q >x ) = (3) 

and "constraint" equations for A 

X,t = t[P 2 + P 2 x + e 2P (Q 2 t + Q 2 x )} (4) 
X, x = 2t(P >x P, t + e 2P Q >x Q it ) (5) 

(here >a = 8/ da). The constraint equations determine A once P and Q are known. 
The integrability conditions for the constraint equations are satisfied as a consequence 
of the evolution equations. Since the evolution equations do not depend on A there is 
essentially a complete decoupling of constraints from evolution equations. Therefore, 
for the purposes of this paper we will treat only equations (j2HSI)- The only restriction 
that the constraints place on initial data for equations (J2HH1) is the following: since A at 
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x = is the same as A at x = 2ir } it follows that the integral from to 2ir of the right 
hand side of equation (JSJ must vanish. We require that this restriction is satisfied by 
the initial data for equations (J2HHI) and then these equations insure that the restriction 
is also satisfied at subsequent times. 

The singularity is at t — 0. It is often helpful to introduce the coordinate r = — Int. 
Thus the singularity is approached as r — > oo. In terms of this coordinate, the evolution 
equations (j2HSl) become 

P,tt ~ e 2P Q 2 T - e~^P, xx + e 2 ^Q 2 x = (6) 
Q, rr + 2P T Q tT - e~ 2r {Q, xx + 2P tX Q tX ) = (7) 

Equations (jHEJ) are the equations of motion corresponding to the Hamiltonian 
H = Hk + Hy where 

H K = irdx(7rl + e- 2P ir 2 Q ), (8) 



2 Jo 

H v = \e-^ rdx(P? x + e 2P Q%). (9) 
2 Jo 

One can also put the evolution equations (J2HHJ) in characteristic form: define null 

coordinates £ = t + x and r\ = t — x and characteristic variables A, B, C and D given by 

a = a + n)p, v , (io) 

B=(£ + 77)F € , (11) 

C=(£ + v)e P Q,r ] , (12) 

D = {i + r 1 )e p Q ti . (13) 

Then equations (J21ISI) become 

A,s = (Z + V r 1 [CD + (A-B)/2] (14) 

S^^ + ^-^CD+^-A)^] (15) 

C,t = (C + r^)- 1 [—AD + (C - D)/2] (16) 

A, = (e + vr'l-BC + (D — C)/2] (17) 

As we will see in the next section, each of these different ways of writing the 
evolution equations has its uses for numerical simulations. 



3. Numerical methods and results 



We now consider numerical simulations of the evolution equations for P and Q. Since 
we wish to examine the behavior as the singularity is approached, it is most convenient 
to simulate equations (jEHZJ) and examine the behavior for large r. One simple numerical 
method is to use centered differences for spatial derivatives and the iterated Crank- 
Nicholson (ICN) method^H] for time evolution. (Spikes were first found0 using the 
symplectic method to be described below. However, the ICN method has been used[TZ| 
to simulate Gowdy spacetimes with spatial topology S 2 x S 1 ). The ICN method is in 
the class of so called "predictor-corrector" methods: first taking a crude approximation 
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Figure 1. Plot of P vs x for vq = 5 at r = 10. Here < x < tt 



to the evolution and then refining it. The symplectic method is in the class of "operator 
splitting" methods: splitting the equations of motion into two pieces, each of which is 
evolved separately. Here centered differences means that for a quantity F at spatial grid 
point i we have F x -> (F i+1 - F i . 1 )/(2Ax) and F )XX -> (F i+1 + F^ x - 2F t )/(Axf. The 
ICN method is the following: we put the equations in the form S jT = W (i.e. first order 
in time) by making P T and Q jT extra variables. Then we evolve the variables S from 
time step n to time step n + 1 using 



We choose S n as our initial guess for S n+ and use this guess in the right hand side of 
equation (fTHj) to produce a better guess. This process is iterated 3 times to evolve to 
the next time step. 

We wish to consider small scale structure generated by the approach to the 
singularity. For that reason we choose initial data that itself has as little small 
scale structure as possible. Following references^ El we use the data P = 0, P r = 
■Uocosx, Q = cosx,Q jT = at time r = 0. (Here Vq is a constant). As argued in 
reference jH] the type of fine structure generated by this family of initial data should be 
generic. Figures ((TJ and (J2J) show the evolution to r = 10 of these data for vq = 5. 

From these plots it is apparent that there are narrow upward pointing features in 
P at x ~ 0.36, x ~ 1.2 and x ~ 2.9. In addition, there is a narrow downward pointing 
feature in P at x ~ 1.66 which coincides with the location of a narrow feature in Q. For 
reasons that will become apparent in section 5, we adopt the notation of reference [IT] 
and call the upward pointing features in P "true spikes" and the downward pointing 
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Figure 4. Plot of Q vs x for vq — 5 at r — 10. This is a closer look at the true spike 
at x w 0.36 





Figure 7. Plot of P vs x for vq = 5 at r = 10 (solid line) and r = 15 (dashed line). 
This is a close look at the true spike at two different times. Note that the spike is 
higher and more narrow at the later time. 
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Figure 8. Plot of Q vs x for vq = 5 at r = 10 (solid line) and r = 15 (dashed line). 
This is a close look at the false spike at two different times. Note that the spike is 
higher and more narrow at the later time. 

features in P along with narrow feature in Q a "false spike." 

Figures Q and (JH) provide a closer look at the true spike at x « 0.36 while figures 
© and (jnj) do the same for the false spike. Note that Q has an extremum at the position 
of the true spike. 

The spikes become steeper and narrower as the singularity is approached. Figure 
(7) shows P at the true spike at x ~ 0.36 at r = 10 and r = 15. Figure (jHJ shows Q at 
the false spike at r = 10 and r = 15. 

The ICN method that we have described is second order accurate. However, one 
might want higher order accuracy, and this is provided by the symplectic methods 
used injS]. The symplectic method^Hj works as follows: begin with a system whose 
equations of motion come from a Hamiltonian of the form H = Hk + Hy where the 
sub-hamiltonians Hk and Hy each correspond to a system that can be solved in closed 
form. (Note that this condition is satisfied by the Gowdy evolution equations (0IZ|) 
where the corresponding sub-hamiltonians are given in equations (JEKI))- Let Wa'(At) 
be the evolution operator that corresponds to evolving the system with Hamiltonian Hk 
for a time At. Correspondingly define Uy(Ar). Now consider the operator W( 2 )(At) 
defined by 

W( 2 )(Ar) = U K (Ar/2)Uy(Ar)U K (Ar/2). (19) 

In words, U^){Ar) evolves using Hk for half a time step, then using Hy for a full time 
step and then using Hk for half a time step. One can show[T£] that U{2) is a second 
order accurate operator for the full system described by H = H K + Hy. Furthermore, 
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using U12) one can construct evolution operators of any accuracy one wishes. The fourth 
order accurate operator Uu) is given by 



from lower order W( n ) in similar ways. 

To take advantage of nth order accuracy in time, the codes must also have nth 
order accuracy in space. Thus the second order accurate centered differences must be 
replaced by more complicated finite differences that are nth order accurate. There is 
another advantage besides higher order accuracy in using the symplectic method to 
study the approach to the singularity: as the singularity is approached, Hy appears to 
become small compared to Hk- If Hy were completely negligible, then the operator W( 2 ) 
(and therefore and all higher order W( n )) would be exact. Thus as the singularity is 
approached, the symplectic method may be even more accurate than one might guess 
from the formal order of accuracy. (This fact has been used in simulations of Mixmaster 
spacetimes [TTJ]. ) 

A numerical method will be accurate only if there are enough spatial points to 
resolve all features. However, the spikes grow ever narrower; thus for any fixed resolution 
there will come a time when that resolution is not sufficient to resolve the spikes. Indeed 
the simulations of E] are often run so long that some spikes are not resolved; however 
it is argued in[Hj that this does not affect the accuracy of the simulation in regions away 
from the spikes. 

To resolve the spikes for a long time, it is better to use a grid that does not have a 
fixed resolution. One way to do this is to use the technique of adaptive mesh refinement 
(AMR).[20J This technique adds extra mesh points in places where the existing mesh is 
not sufficiently fine to resolve the features. AMR was applied to the Gowdy spacetimes 
by Hern and Stewart. |2"T] Unfortunately, AMR codes are quite involved and difficult to 
write. Moreover, the full machinery of AMR may not be needed since it appears that 
in a given Gowdy spacetime there is a limited number of spikes and that their positions 
change very little as the singularity is approached. One can then construct a mesh that 
is finer in the places where one knows the spikes will be.|22j 

Alternatively, one can study a single spike using a characteristic method. 23 a Here, 
the initial data surface is an outgoing light ray and one evolves along ingoing light rays. 
One can choose the last grid point to correspond to the light ray that hits the center 
of the spike at the singularity. As a consequence of this choice, the grid shrinks as the 
evolution proceeds in such a way that the grid is always of an appropriate size to resolve 
the spike. 

4. Analytical approximations 

We now turn to analytical approximations to the equations of motion (jBHZ|). Our 
treatment closely follows that of[6|. One obvious approximation is the following: since 



where s = 1/(2 
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the approach to the singularity is r — > oo, it is natural to neglect terms in equations 
(00 proportional to e~ 2r . Dropping such terms results in the following equations: 

P,rr - e 2P Q 2 T = 0, (21) 
Q, rr + 2P )T Q jT = (22) 

These equations are called the velocity term dominated (VTD) equations. (Note that 
the VTD equations are the equations of motion corresponding to the sub-hamiltonian 
H K ). A solution of equations (00 is called asymptotically velocity term dominated 
(AVTD) if there is a solution of the VTD equations that it approaches as r — > oo. 

The VTD equations can be solved in closed form with the general solution given 

by 

P — p + lnfcoshwr + cos-^sinhur], (23) 
^ e~ p sin ib tanh vt . iS 

Q ^ + l + co^ t anh OT (24) 
Here i> > and the quantities p, q, v and ib are functions of x. The behavior of the false 
spikes follows from equations (|23M24j) . For costp 7^ — 1 the behavior of these solutions 
as t —> oo is P —> p + vt and Q —> q where p and q are functions of x given by 
p = p + ln[(l + cos^)/2] and q = q + e~ p sin^/(l + cos^)- Let x\ be a point where 
cosV 1 = — 1. Then clearly there is some steep behavior near x\. Let a subscript 1 denote 
the value of the function at x% and let ib\ denote ib >x at x\. Then near x% and for large 
r it follows from equations (|23l - 12*^) that 

P w pi - uir + ln[l + (^i') 2 e 2,,ir (x - xi) 2 /4] (25) 

v ~ Qi 5 — (26) 

2e- 2 ^ + (^ l ') 2 (x-x 1 ) 2 /2 



Equations (j25H26j) provide an analytic approximation to the behavior of false spikes. 
Figures ([51110)) show the behavior of P and Q respectively for these formulas. Here 
Pi — 5, qi — 0, v\ = 0.5, x± = 0, ipi = 1 and r = 10. 

To treat the behavior of the true spikes, we first consider the validity of the VTD 
approximation by using the "method of consistent potentials" (MCP).j^ That is, we 
assume the large r behavior P — > p + vt, Q —>■ q and see whether the terms in equations 
(0-0 that we have neglected are in fact negligible. The term e 2 ( p ~ T >Q 2 x has behavior 
e 2[p+(i;-i)r]^_ Thus this term is negligible only if v < 1. The numerical simulations 
of j3J show that in fact v evolves to a value less than 1 except at the true spikes. To 
find an explanation for this behavior we must consider what happens when v is initially 
greater than 1. We are still justified in using equation which yields Q T = 7iQe~ 2P 
where tiq is a function of x. This turns the term —e 2P Q 2 T in equation (0 into — e~ 2P iiQ 
which can be neglected since P is growing. We can also approximate Q{r,x) by q{x). 
Thus we are led to approximate equation (0 by 

P TT + e* p - T % = 0. (27) 

This equation can be solved in closed form to yield 

P = p + t — In [cosh wt — cos (f) sinh wt] . (28) 
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Figure 9. Plot of P vs x for the analytic approximation, equation Ij25(l for a false 
spike 
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Figure 11. Plot of P vs x for the analytic approximation, equation (|29(l for a true 
spike 

Here w > and p, w and </> are functions of x and we have q yX = e~ p w sin <p. For cos <p ^ 1 
the large r behavior of this solution is P —> p + vt where p = p — ln[(l — cos0)/2] and 
v = 1 — w. Let x 2 be a point where cos0 = 1 and use a subscript 2 to denote the value 
of the function at x 2 . Then near x 2 and for large r if follows from equation (}28|) that 

P^p 2 + (l + w 2 )r - ln[l + {<p 2 ' fe 2w2T (x - x 2 ) 2 /4] (29) 

Equation (|29|) provides an analytic approximation to the behavior of true spikes. Figure 
(111}) shows the behavior of P for this formula. Here p 2 = 0, w 2 = 0.5, x 2 = 0, 2 ' — 1 
and r = 10. 

The dynamics of spike formation can be understood as follows: Define the 
"potentials" V\ = TCQe~ 2P and V 2 = e, 2( - p ~ T ^Q 2 x . Then the behavior of P can be regarded 
as the the dynamics of a particle in these potentials. If P T < then a "bounce" off 
of V\ will give rise to a transition P T — > — P T which will make P T positive. If P T > 1 
then a bounce off of V 2 will give rise to a transition P r — > 1 — P T and thus P T < 1. 
Eventually, after a finite number of bounces P T will be in the range < P T < 1. This 
occurs provided that both potentials Vi and V 2 are nonvanishing. At a point xi where 
Vi = a false spike can form since P T < is allowed. At a point x 2 where V 2 = a 
true spike can form since P r > 1 is allowed. 

We now consider the validity of the approximation (j2*9*|) . Using the MCP we find 
that the term e~ 2r P yXX goes as e 2( - w ~~ 1 ^ T . This term is thus negligible provided that w < 1 
or equivalently that P T at the center of the spike is less than 2. This leaves open the 
question of what happens to "high velocity" spikes, i.e. those for which P T at the 
center of the spike is initially larger than 2. This question was answered recently in the 
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numerical work of (22]. The simulations of (22] show that initially high velocity spikes 
are forced by the term e~ 2r P^ xx to lower velocity. In the end, either the spike itself 
disappears or the velocity at the center of the spike is brought to the range 1 < P r < 2. 

5. Mathematical results 

Here we consider two mathematical questions about spikes: (1) are spikes geometrical 
features or coordinate artifacts? and (2) to what extent are solutions with spikes 
rigorously proven to exist? This section is essentially a summary of the results of [llj . 
The reader should see that work and references therein for a more detailed treatment. 

So far we have presented the spikes in terms of the behavior of P and Q as functions 
of r and x. That is, we have considered the coordinate dependence of metric components. 
One might then wonder whether the spikes are simply artifacts of the coordinate system 
in which the metric is expressed. However, the Gowdy coordinates and metric 
components themselves have geometric meaning in terms of the symmetries of the 
spacetime. The coordinate t is defined by the area of each symmetry T 2 being 4n 2 t. The 
coordinate x is the harmonic function conjugate to t (see Chapter 3, problem 2 of |24|). 
The vector fields (d/dy) a and (d/dz) a are Killing fields. Furthermore, P and Q can be 
expressed in terms of inner products of Killing fields. Nonetheless, there is a restricted 
coordinate freedom that preserves this relation between coordinates and geometry. In 
particular, consider the transformation that simply switches the coordinates y and z. 
This leaves the metric in the form of ([TJ but changes the pair (P, Q) to (P, Q) given by 

e ""=Q^' (30) 
Q_ 

Q 2 + e~ 

Following ^1] we will use the notation (P, Q) = I(P, Q) where the "I" is for inversion. 
Now suppose that (P, Q) is a solution of the VTD equations without spikes and which 
therefore has the large r behavior P — > p + vr and Q — > q. Then Q has the large r 
behavior 

<2 -> q2 + e q 2p e -2 VT - ( 32 ) 

Therefore Q has a false spike at points where q vanishes. Thus false spikes are not 
geometric quantities since they can be made to appear and disappear by a coordinate 
transformation. 

In contrast true spikes can be shown to be true geometric objects. This is done 
by examining the behavior of curvature invariants at the spikes as the singularity is 
approached. The asymptotic behavior of curvature invariants as the singularity is 
approached is different at a true spike than at nearby points. 

We now turn to a different type of transformation that takes a solution (P, Q) of 
equations (jMZj) and produces a solution (P, Q) that represents a different spacetime. 



Q = T^~^- (31) 
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The transformation is given by 

P = r-P, (33) 

Q,t = -e 2iP ~ T) Q, x , (34) 

Q,x = ~e 2P Q, T . (35) 

This transformation is essentially the analog for Gowdy spacetimes of the Kramer- 

Neugebauer transformation [2*5] for stationary axisymmetric spacetimes. The 



integrability conditions for equations (|34H35[) are satisfied as a result of the equations 
of motion (jBHZj). Following 11] we will use the notation (P,Q) = GE(P, Q). Here the 
"GE" stands for "Gowdy to Ernst." For our purposes what is important is that if (P, Q) 
has a false spike then GE(P, Q) will have a true spike. Thus from a solution (P, Q) with 
no spikes at all we can generate a solution GE(I(P, Q)) with true spikes. The problem 
of proving existence of solutions with spikes has thus been reduced to the problem of 
proving existence of solutions with no spikes. 

How then does one prove existence of solutions with no spikes? Here we simply 
sketch what is done and refer the reader to (HI Ej for details. At first it might seem that 
what is needed is a global existence theorem. However, spikes are essentially asymptotic 
behavior as the singularity is approached. Thus a local existence theorem will suffice 
provided that it is local in a neighborhood of the singularity. What is needed for such 
a theorem is to put the equations in Fuchsian form 
du -* 

t— +N(x)u = tf(t,x,u,u tX ) (36) 

for a system described by the vector valued function u. One needs certain properties 
for the function / and matrix N. Using these properties, one proves that there exists 
a solution in a neighborhood of t = and with u — at t — 0. To write the Gowdy 
equations in Fuchsian form, one must first correctly guess (based on equations (|23H24p ) 
functions that have the t — > behavior of P and Q. One then writes P and Q as these 
functions plus remainder terms. One then rewrites the Gowdy evolution equations (J2K3J) 
as a Fuchsian system for the remainder terms. In this way existence of P and Q with 
the appropriate asymptotic behavior is proved. The P and Q obtained in this way do 
not have spikes. However, applying the transformations I and GE to this P and Q one 
obtains a rigorous proof of existence of Gowdy spacetimes with spikes. 

6. Conclusions 

When spikes were first found by Berger and Moncrief [5] they seemed quite mysterious. 
Now, about a decade later, we have a very good understanding of this phenomenon. 
A variety of numerical simulations have allowed us a close look at the process of spike 
formation and at the properties of the spikes. Analytical approximations give us an 
understanding of how the spikes form as well as an approximate description of their 
shape and time development. Mathematical results have shown that the false spikes 



The fine structure of Gowdy spacetimes 



15 



are mere coordinate artifacts while the true spikes are real geometric phenomena. In 
addition it has been possible to prove rigorously that solutions with spikes exist. 
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